function [char_fiber_length, char_pore_size] = CollagenAutocorrelation(image_double, pixel_size, resizing_factor, weiner_filt, filter_size, top_hat_filt, strel_disk, shift, BW_threshold, center_zero)

%% COLLAGEN AUTOCORRELATION ANALYSIS
% Authors: Francois Bordeleau, Matthew Zanotelli, Paul Taufalele
%   last edit 2018.7.13
% Purpose:  take advantage of autocorrelation analysis to estimate the
%           characteristic feature length and pore size of collagen confocal
%           reflectance images.  
% Inputs:
%       image_double - a square matrix of type double - stores the image to
%       be analyzed
%       pixel_size - scalar - # of pixels per micron
%       resizing_factor - scalar - to shrink image in order to speed up
%       analysis.  I have found 0.90 for 1024x1024 images to work well.
%       weiner_filt - scalar - set to 1.0 to use the weiner filter 
%           filter_size - scalar - pixel size for weiner filter
%       top_hat_filt - scalar - set to 1.0 to use top hat filter 
%           strel_disk - scalar - pixel size for strel disk used in top hat
%           filter
%       shift - scalar - size of shift used for autocorrelation (at least
%       2x the size of desired feature
%       BW_threshold - scalar - intensity threshold for generating binary
%       image.  I have found 0.01 to work well, but will need to be
%       adjusted each time
%       center_zero - scalar - set to 1.0 to set the center autocorrelation
%       pixels to 0.  This helps remove artifacts generated from the
%       fourier transform.
% Outputs: 
%       char_fiber_length - scalar - characteristic feature length in um
%       generated from autocorrelation analysis of weiner/tophat filtered
%       image
%       char_pore_size - scalar - characteristic pore size in um generated
%       from autocorrelation analysis of BW image
%       GRAPHS - displays the autocorrelation images of the filtered/BW
%       images with the contour used to calculate the char_fiber_length and
%       char_pore_size.  Also displays a graph 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input calculation parameters
% Microsocpe calibration (in um/pixel); based on microscope and objective
% used for imaging
calibration = (1/pixel_size);
% Resizing factor
% This is used to increase the speed of analysis; i.e. if you use 0.5, a
% 1024x1024 image will become 512x512
s_Im = resizing_factor;

%% Image load sequence
image = image_double;

%% Applying desired image filtering steps with user input values
if weiner_filt == 1
    image = wiener2(image,[filter_size filter_size]);
end

if top_hat_filt == 1
    tophatFiltered = imtophat(image,strel('disk',strel_disk));
    image = imadjust(tophatFiltered);
end

image = imresize(image,s_Im); % resize the image

% Sum of all pixel intensities in 'image'
total1=sum(sum(image)); 
% Normalize the autocorrelation.
total_poresize=sum(sum(max(image(:))-image));                        

image1 = max(image(:))-image; % inverts pixel intensities
image1 = image1-min(image1(:)); % subtract min pixel intensity from image
image1 = image1-mean(image1(:)); % subtract mean pixel intensity from image
BW = im2bw(image1,BW_threshold); % generate BW image for pore size autocorrelation

figure % display BW image
imshow(BW);

[N M Rgb]= size(image1);
tau = 1/(M*N); 
acf=zeros([2*shift+1 2*shift+1]); % Result matrice initialization

% autocorrelation
nfft = 2^nextpow2(2*N-1); % some sort of padding for fft
acfFL = fftshift(ifft2(fft2(image1).* conj(fft2(image1))));
acfPS = fftshift(ifft2(fft2(BW).* conj(fft2(BW))));
% Normalize the autocorellation to min in image
AcorrFL=acfFL(N/2-shift:N/2+shift,M/2-shift:M/2+shift)-min(min(acfFL(N/2-shift:N/2+shift,M/2-shift:M/2+shift))); 
Acorr_poresize=acfPS(N/2-shift:N/2+shift,M/2-shift:M/2+shift)-min(min(acfPS(N/2-shift:N/2+shift,M/2-shift:M/2+shift))); 
 
AcorrFL=AcorrFL/max(AcorrFL(:))*255;    % Then setting a 8bit pseudo scale
Acorr_poresize=Acorr_poresize/max(Acorr_poresize(:))*255;
maxA=max(AcorrFL(:));
maxB=max(Acorr_poresize(:));
if center_zero == 1
    AcorrFL(AcorrFL==maxA)=0;               % The fundamental frequency can be removed since it
    Acorr_poresize(Acorr_poresize==maxB)=0; % does not contain information about the features
    maxA=max(AcorrFL(:));
    maxB=max(Acorr_poresize(:));
end

eA=(1/exp(2))*(maxA);   % using the 1/e^2 criteria to find the feature edge.
eB=1/exp(2)*(maxB);   % this criterria should contain 99% of the frequency 
                     % describing the lenght scale of the image

                    
%% Image generation

% initialization required for contour function
x = 1:101;
y = 1:101;
[X,Y] = meshgrid(x,y);
% full scale image
Ls = figure;
imagesc(AcorrFL,[0 maxA]);figure(Ls);
colormap(gray)
hold 'on';
[Ca] = contour(X,Y,AcorrFL,[eA,eA],'red');     % finding the 1/e^2 contour
[Cb] = contour(X,Y,AcorrFL,[maxA/2, maxA/2],'red'); % finding the I(0)/2 contour

Lp = figure;
imagesc(Acorr_poresize,[0 maxA]);figure(Lp);
colormap(gray)
hold 'on';
[Da] = contour(X,Y,Acorr_poresize,[eB,eB],'red');      % finding the 1/e^2 contour
[Db] = contour(X,Y,Acorr_poresize,[maxB/2, maxB/2],'red');  % finding the I(0)/2 contour

% computing the average radius of the area defined by the 1/e^2 limit
% Ma = mean(sqrt((Ca(1,2:Ca(2,1))-shift-1).^2+(Ca(2,2:Ca(2,1))-shift-1).^2))*calibration/s_Im*2;
Ma = 2*sqrt(polyarea(Ca(1,2:end), Ca(2,2:end))/pi)*calibration/s_Im*2;
% Mb = mean(sqrt((Da(1,2:Da(2,1))-shift-1).^2+(Da(2,2:Da(2,1))-shift-1).^2))*calibration/s_Im*2;
Mb = 2*sqrt(polyarea(Da(1,2:end), Da(2,2:end))/pi)*calibration/s_Im*2;  % Diameter of circle w/ same area as contour

figure(Ls)
annotation('textbox',[0.6 0.7 0.1 0.1],'String',...
    [texlabel('omega = ') num2str(Ma,4) texlabel(' mu') 'm'],...
    'FontSize',14,'linestyle','none')

figure(Lp)
annotation('textbox',[0.6 0.7 0.1 0.1],'String',...
    [texlabel('omega = ') num2str(Ma,4) texlabel(' mu') 'm'],...
    'FontSize',14,'linestyle','none')

% Numerical values are stored is this structure
% 2D plot
Lp = figure;
plot(AcorrFL(shift+1,:));figure(Lp);

%% Output values
char_fiber_length = Ma;
char_pore_size = Mb;
length = {'Char fiber length:',Ma; 'Char pore size:',Mb};
 
end

 